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Abstract 

The NAIAD experiment (Nal Advanced Detector) for WIMP dark matter 
searches at the Boulby Underground Laboratory (North Yorkshire, UK) ran from 
2000 until 2003. A total of 44.9 kgx years of data collected with 2 encapsulated 
and 4 unencapsulated Nal(Tl) crystals with high light yield were included in 
the analysis. We present final results of this analysis carried out using pulse 
shape discrimination. No signal associated with nuclear recoils from WIMP 
interactions was observed in any run with any crystal. This allowed us to set 
upper limits on the WIMP-nucleon spin-independent and WIMP-proton spin- 
dependent cross-sections. The NAIAD experiment has so far imposed the most 
stringent constraints on the spin-dependent WIMP-proton cross-section. 
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1. Introduction 

The UK Dark Matter Collaboration (UKDMC) has been operating various detec- 
tors for non-baryonic dark matter searches at the Boulby Underground Laboratory 
(North Yorkshire, UK) for many years. Limits on the flux of weakly interacting mas- 
sive particles (WIMPs) - the primary candidate for the non-baryonic dark matter that 
may constitute up to 90% of the mass of the Galaxy, were set using data from the first 
Nal(Tl) detector |2j and later improved with an array of several Nal(Tl) crystals 
(NAIAD) 3 . Pulse shape analysis (PSA) was applied to the data to distinguish be- 
tween slow scintillations arising from background electron recoils and fast scintillations 
due to nuclear recoils, which are expected from WIMP-nucleus interactions jl]. 

The advance of the UKDMC Nal(Tl) experiment was blocked for a few years (1997- 
2000) by the discovery of a fast anomalous component in the data from several encap- 
sulated crystals jS]. These events were faster than typical electron recoil pulses and 
faster even than nuclear recoil pulses jSlEl- Similar events were also seen by the Saclay 
group [7j . These events were later shown to be due to implanted surface contamination 
of the crystal by an alpha-emitting isotope from radon decay |HJ E] • 

NAIAD was constructed in 2000-2001 with 5 unencapsulated crystals to allow better 
control of the crystal surface, to reduce the background due to surface events and to 
improve the light collection In 2002 two more encapsulated crystals were added 

to the array. 

Despite the limited sensitivity of inorganic scintillators to WIMP interactions, in- 
terest in the Nal(Tl) crystals has remained high for the past several years because of 
the DAMA group's claim of an annual modulation in their Nal(Tl) array consistent 
with the expected signal from WIMP-nucleus interactions with a specific set of WIMP 
parameters [TT] . 

Nal has the advantage of having two target nuclei with high and low masses, thus re- 
ducing uncertainties related to nuclear physics calculations. The detectors are sensitive 
to both spin-independent and spin-dependent interactions. The NAIAD experiment is 
complementary to other dark matter experiments at Boulby, such as ZEPLIN ^2] (liq- 
uid xenon and two phase xenon detectors) and DRIFT (time projection chamber 
with directional sensitivity). The array of Nal(Tl) detectors was also used as a di- 
agnostic array to study background and systematic effects for other experiments at 
Boulby. 

First results from the NAIAD experiment (10.6 kgxyears exposure from 2000-2001) 
were published in Ref . j3] . Preliminary limits from an extended set of data (combined 
with the previous set) with total exposure of 20.3 kgxyears (2000-2002) were set in 
Ref. jHj. In this paper we present results from the NAIAD experiment based on all 
data sets collected in 2000-2003. 

2. The NAIAD experiment 

The NAIAD array was sited in the underground laboratory at Boulby mine at a 
vertical depth of 1070 metres or 2805 m w. e. [T5] . It consisted of 7 Nal(Tl) crystals 
from different manufacturers (Bicron, Hilger, VIMS and Crismatec-Saint-Gobain) with 
a total mass of about 55 kg. Two detectors contained encapsulated crystals, while 5 
other crystals were unencapsulated. To avoid their degradation by humidity in the 
atmosphere, the unencapsulated crystals were sealed in copper boxes filled with dry 
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nitrogen. 

One of the crystals (DM70-Saclay) of similar design to those used by the DAMA 
group in its 100 kg array, was previously running at Modane and then moved to Boulby 
in 2001. When encapsulated, it showed the anomalous fast population of events 
which were explained as being due to implanted surface contamination of the crystal 
by an alpha-emitting isotope from radon decay jHUHl- It was de-encapsulated in 2002, 
polished and was running as part of the NAIAD array in 2002-2003. No fast population 
was seen in the crystal after de-encapsulation [TBI ITTj. 

Each crystal was mounted in a 10 mm thick solid PTFE (polytetraffuoroethylene) 
reflector cage and was coupled to light guides. The two 4-5 cm long quartz light 
guides were also mounted in the PTFE cages and were coupled to 5 inch diameter 
low background photomultiplier tubes (PMTs), ETL type 9390UKB. Only selected 
low background materials were used in the detector design including oxygen-free high- 
conductivity copper and PTFE. 

Temperature control of the system was achieved through copper coils outside the 
copper box. Chilled water was constantly pumped through the coils maintaining the 
temperature of the crystals stable to within a fraction of a degree during a single 
run. The temperature of the crystal, ambient air, water in the pipes and copper 
was measured by thermocouples. If, for any reason (for example, chiller failure), the 
variation of crystal temperature exceeded the predefined limit, the data from these 
periods were not included in the analysis. If any changes in the experimental set-up 
resulted in a change of the mean temperature of the crystal, the data from different 
runs were not combined together, but assumed to come from different experiments, so 
only the resulting limits were combined. 

Pulses from individual PMTs were integrated using a buffer circuit and then digi- 
tised using an Acqiris CompactPCI based data acquisition (DAQ) system. Note that 
previous results jS] were obtained with a digital oscilloscope and Labview-based DAQ 
software. The digitised pulse shapes (5 fis total digitisation time, 10 ns digitisation 
accuracy) were passed to a computer running Linux OS and stored on disk. The gain 
of the PMTs was set to give about 2.5 mV per photoelectron. Low threshold discrim- 
inators were set to about 10 mV threshold, which corresponded to about 1-1.5 keV of 
electron equivalent energy. A software threshold was set at 4 keV since no pulse shape 
discrimination was observed below this energy [3]. 

Copper boxes containing the crystals were installed in lead and copper "castles" , to 
shield the detectors from background due to natural radioactivity in the surrounding 
rock. Wax and polypropylene neutron shielding (about 10 g/cm 2 ) around the castles 
was installed in Spring 2002 and all data presented here (2002-2003) were collected 
with neutron shielding in place. 

The light yield of the crystals was obtained from measurement of the single pho- 
toelectron pulses and the energy calibration with a 122 keV 57 Co gamma-ray source. 
The crystals had light yields of 4.6-8.4 pe/keV. The large difference in the crystal light 
yield is not surprising taking into account the different origin (manufacturers) of the 
crystals, encapsulation, reflecting material, size of the light guides and PMT response. 
The light yield was checked every 2-4 weeks and was found to be stable within 10% for 
all crystals over long running periods (several months and more). The longest operated 
unencapsulated crystal, running since the summer of 2000, showed a degradation in 
light yield of no more than 10%. Small variations in the light yield do not affect the 
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energy threshold for the data analysis since it is significantly higher (4 keV) than the 
effective hardware threshold (typically less than 2 keV). Possible short term variations 
in the light yield, as well as any electronic problems (stability of the discriminator 
thresholds, high- voltage power supplies etc.) were monitored by plotting the energy 
spectra of gamma calibration data and real data for each day of running and remov- 
ing suspicious data from the analysis. Any variations in the experimental conditions, 
mentioned above, would manifest themselves in the energy spectra. 

3. Analysis procedure 

Final analysis was performed on the sum of the pulses from the two PMTs attached 
to each crystal. The parameters of the pulses from each PMT were used to apply 
various cuts to remove the noise events from the analysis. Our standard procedure of 
data analysis involved fitting a single exponential to each integrated pulse in order to 
obtain the index of the exponential, r, the amplitude of the pulse, A, and the start 
time, t s . The scintillation pulses from nuclear and electron recoils were shown to be 
nearly exponential in shape jl]. In addition to these, the mean time of the pulse (mean 
photoelectron arrival time), as an estimate of the time constant, x 2 °f t ne fit, number 
of photoelectrons and the energy were also calculated for each pulse. Conversion of the 
pulse amplitude to the number of photoelectrons and the energy was achieved using 
pre-determined conversion factors, which come from energy and single photoelectron 
calibrations. These calibrations were performed typically every 2-4 weeks to ensure 
that there were no changes in the PMT gains or crystal degradation. Calibration with 
various gamma-ray sources provided the energy resolution of crystals [3 . 

For each run (or set of runs) the "energy - time constant" (E — r) distribution 
was constructed. If all operational settings (including temperature) were the same for 
several runs, the (E — r) distributions for these runs were summed together. Data 
collected at crystal temperatures outside a predefined range were removed from the 
analysis. Each event stored on the computer disk, had a temperature record in it. 
The mean daily temperature of each crystal was plotted as a function of time, and 
the mean temperature for a particular data set was calculated. The most common 
reason for a temperature to fluctuate significantly (several degrees) from the mean 
value was chiller failure. Data collected without cooling were removed from the analysis. 
When a broken chiller was replaced or repaired, it was impossible to achieve the same 
crystal temperature. In this case we considered the two data sets collected at different 
temperatures as independent and did not combine the data but only the final limits 


To reduce PMT noise and, particularly, events where a spark (flash) in the dynode 
structure of one PMT was seen by both PMTs, various cuts described in details in 
[TB| ITS] were applied. These cuts are an improvement upon those used in Ref. [3], 
having been shown to reject spurious events more efficiently while retaining a larger 
fraction of genuine scintillation pulses [TBI EE]- The noise cuts were developed based 
on the experiment with two PMTs mounted face to face without the crystal between 
them. They were shown to cut similar fraction of scintillation pulses from electron and 
nuclear recoils (to within 2%) [TBI ITS]. The cut efficiency for nuclear recoil events from 
the surface calibration ranges from 84% at 4-6 keV up to almost 100% above 10 keV 
[TBI ITS] for one of the crystals. As these cuts were found to be similar (to within 2%) 
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for electron and nuclear recoils, the cut efficiencies for each particular crystal and data 
set were determined from the daily calibrations with the 60 Co gamma source. This also 
allowed the monitoring of the stability of these cuts. 

For any small energy bin (1 keV width, for example), the time constant distribution 
was approximated by a Gaussian in ln(r) (log(Gauss) function) [21 El El UHj with three 
free parameters: mean time constant r a , width w and normalisation factor N a . In 
experiments where a second population is seen (for example, nuclear recoils from a 
neutron source or possible WIMP-nucleus interactions), the resulting r-distribution 
was fitted with two log(Gauss) functions with the same width w, since the width is 
determined mainly by the number of collected photoelectrons. 

The aim of the analysis procedure was to find a second population of events in the 
r-distributions or to set an upper limit on its rate. To reach this, the r-distributions 
for all energies of interest (2-40 keV) were obtained for all crystals with gamma-ray 
( 60 Co) and neutron ( 252 Cf) sources. Photons from high-energy gamma-ray sources 
produce Compton electrons in the crystal volume similar to those initiated by gamma- 
ray background in the crystal. Neutrons collide elastically with the nuclei of the crystal 
giving nuclear recoils similar to those expected from WIMP-nucleus elastic scattering. 
Calibrations were performed on all crystals at the surface prior to moving the crystals 
underground and were described in detail in Refs. E] ■ 

The crystal temperature and the light yield are the most important parameters 
which affect the time constant distributions and, finally, the results of the experiment. 
The temperature affects the pulse shape - mean time of the pulse and the index of an 
exponential. Increasing temperature causes the mean time (and the index of the fitted 
exponential) to decrease at a rate of about 5 ns per 1°C. If the temperature of the 
crystal is not stable, then combining data collected over a long period of time, we will 
broaden the time constant distribution. The same effect is observed if the light yield is 
not stable. To improve the quality of the data, the crystal temperature was controlled 
and monitored. To compensate for any small variations of the crystal temperature and 
light yield over long periods of time (several months) the gamma calibration of the 
crystal was done every day for 2 hours with the 60 Co gamma-ray source. In this way, 
any variation in the mean time constant and width of the time constant distributions 
of the real data, due to the effects mentioned above, is repeated in similar distributions 
constructed with gamma calibration data and used to fix some parameters in the fitting 
procedure (see Ref. |3] for a full description of the data analysis procedure). Note that 
the ratio of mean time constant for nuclear recoils to that of electron recoils does not 
depend on temperature |3J and was taken from the surface calibrations. 

4. Results and discussion 

The data from 6 crystals were used to set the limits on WIMP-nucleus cross-section 
reported here. One crystal was excluded from the data analysis due to its small mass 
(4 kg) and high background rate (about 15 events/kg/day/keV). Table shows the 
main characteristics and statistics for all detectors for years 2002-2003. Statistics for 
the 10.6 kgxyears of data collected in the period 2000-2001 are to be found in Ref. [3]. 

An energy range from 4 to 10 keV was used in the data analysis. Previous results |3J 
were obtained for the energy range 4-30 keV. The smaller energy range in the present 
analysis was due to the smaller amplitude range set on the 8-bit waveform digitisers. A 
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total of 256 digitisation points for a smaller amplitude range allowed a more accurate 
calculation of the amplitude and time constant for single pulses at low energies and 
better noise suppression. 

Figure n shows typical time constant distributions at 5-6 keV from data (before and 
after cuts) and a calibration run with gamma-ray source (after cuts). PMT noise events 
are seen at small values of time constant before cuts but are absent after the cuts were 
applied. Figure [T] does not reveal any visible difference between data and calibration 
runs in terms of time constant distributions. Limits on the nuclear recoil rate at any 
particular energy were obtained by fitting the measured time constant distribution 
with two log(Gauss) functions having known parameters: mean time constants and 
widths, known from Compton and neutron calibrations. Free parameters were the 
total numbers of electron and nuclear recoils (see Ref. [3] for further details). When 
both positive and negative values for the numbers of nuclear and electron recoils were 
allowed, the best fit values for nuclear recoil rates were either positive or negative 
but normally within 1.5 standard deviations from zero. Total event rates and nuclear 
recoil rates extracted from the best two log-Gaussian fits to the data are shown in 
Figure 121 for two runs with crystal DM74. The error bars for nuclear recoil rates 
correspond to 90% CL. As can be seen, in none of the energy bins does the nuclear 
recoil rate deviate significantly from zero. Hence no contribution from WIMP-nucleus 
interactions were observed in these data. This was found to be true for all crystals. 
Since negative values are non-physical, to set upper limits on nuclear recoil rates we 
restricted the rates to non-negative values and carried out a two component fitting to 
the time constant distributions. To derive an upper limit at 90% CL. we set the value 
of A(% 2 ) = xip ~ Xmin = 2-7. Strictly speaking this A(x 2 ) corresponds to the 90% 
CL. only if the number of searched events is equal to 0. For negative values of best 
fit number of nuclear recoils the value for an upper limit obtained in this analysis is 
higher than suggested in Ref. [20] where a unified approach to confidence intervals 
was discussed for the Gaussian-with-boundary problem. For positive values of best fit 
number of nuclear recoils we used the number of recoils at \ 2 as an estimate for an 
upper limit, which again gave us a higher value than would follow from Ref. |20j . Our 
approach was thus conservative and allowed us to avoid complications with different 
approaches suggested by different authors (see Ref. [2UJ an d references therein for 
discussion). In this way, the upper limits on the nuclear recoil rate were obtained for 
each energy bin and for each crystal. 

The limits obtained on the nuclear recoil rate for each energy bin and each crystal 
were converted into limits on the WIMP-nucleon spin-independent and WIMP-proton 
spin-dependent cross-sections following the procedure described by Lewin and Smith 
[2~T] and used previously in Ref. [3]. Expected nuclear recoil spectra from WIMP- 
nucleus interactions were calculated for a model with spherical isothermal dark matter 
halo with parameters: pd m = 0.3 GeV cm -3 , v Q = 220 km/s, v esc = 600 km/s and 
VEarth — 232 km/s. For the spin-independent case the form factors were computed 
using Fermi nuclear density distribution with the parameters to fit muon scattering 
data as described in Ref. [21]. The WIMP-nucleus scattering cross-section was taken 
proportional to A 2 . For the spin-dependent case a pure higgsino was assumed. The 
spin factors and form factors were computed for sodium and iodine nuclei on the basis 
of nuclear shell model calculations [22] with Bonn A potential. The quenching factors 
(scintillation efficiencies for nuclear recoils with respect to that for electron recoils) 
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were taken as 0.275 for sodium and 0.086 for iodine recoils [3]. 

The limits on the cross-section for various energy bins, targets (sodium and iodine) 
and crystals were combined following the procedure described in Ref. |2lJ and used 
previously in Ref. [3] , assuming the measurements for different energy bins and different 
crystals were statistically independent. 

Figured (Gfc>) shows the NAIAD limits on WIMP-nucleon spin-independent (WIMP- 
proton spin-dependent) cross-sections as functions of WIMP mass based on the data 
described in Table H Also shown in Figure is the region of parameter space favoured 
by the DAMA positive annual modulation signal [TT] . 

Model-independent limits on spin-dependent WIMP-proton and WIMP-neutron 
cross-sections, calculated following the procedure described in Ref. are presented 
in Figure HJ The NAIAD experiment imposes the most stringent constraints so far on 
the spin- dependent WIMP-proton cross-section (see Refs. (231 EH! for a compilation of 
recent results). 

5. Conclusions 

Results from the NAIAD experiment for WIMP dark matter search at Boulby 
mine were presented. Pulse shape analysis was used to discriminate between nuclear 
recoils, possibly caused by WIMP interactions, and electron recoils due to gamma- 
ray background. We obtained upper limits on the WIMP-nucleon spin-independent 
and WIMP-proton spin-dependent cross-sections based on the data accumulated by 6 
modules (44.9 kg x years exposure). The NAIAD experiment imposes the most stringent 
constraints so far on spin-dependent WIMP-proton cross- sect ion. 
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Table 1: Statistics for NAIAD detectors running during years 2002-2003. 



Crystal 


Mass, kg 


Light yield, pe/keV 


Time, days 


Exposure, kg x days 


DM70-Saclay 


10.0 


5.1 ±0.4 


157.2 


1571.7 


DM74 


8.40 


8.4 ±0.4 


298.3 


2505.4 


DM76 


8.32 


4.6 ±0.3 


131.0 


1089.6 


DM77 


8.41 


6.1 ±0.3 


221.9 


1866.2 


DM80 


8.47 


8.0 ±0.5 


312.0 


2642.7 


DM81 


8.47 


6.6 ±0.4 


336.2 


2847.2 


Total exposure 








12523 
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Figure 1: Typical time constant distributions (shown here for crystal DM74) at energies 
5-6 keV: a - for gamma calibration run (after cuts), b - for data before (open circles) 
and after (filled circles) cuts. 
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Measured energy, keV 



Figure 2: Energy spectra after cuts from two runs of one of the crystals (DM74) (filled 
and open circles) and the nuclear recoil rate for various energy bins (filled and open 
squares) with error bars drawn at 90% CL. Energy bins are: 4-5 keV, 5-6 keV, etc. 
The points are shifted with respect to each other along x-axis to avoid overlapping. 
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Figure 3: NAIAD limits (90% C.L.) on WIMP-nucleon spin-independent (a) and 
WIMP-proton spin-dependent (b) cross-sections as functions of WIMP mass: dashed 
curves - previously published limit |HJ, thin solid curves - limit from 2002-2003 data, 
thick solid curves - combined limit (2000-2003). Also shown is the region of parameter 
space favoured by the claimed DAMA positive annual modulation signal (DAMA/Nal- 
1 through DAMA/NaI-4) (closed curve). 
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Figure 4: Model-independent limits (90% C.L.) on spin-dependent WIMP-proton (a) 
and WIMP-neutron (b) cross- sect ions as functions of WIMP mass. The limits were 
derived following the procedure described in Ref. [23J. Dashed curves show the limits 
extracted from interactions with sodium, dotted curves show those from iodine, and 
solid curves show combined limits. 
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